The present invention relates to migration of seismic reflections performed in a high-speed digital computer, and more particularly to suppressing reflection artifacts imposed on a numerical model by artificial boundaries which limit computational size.
For many years seismic exploration for oil and gas has involved the use of a source of seismic energy and its reception by an array of seismic detectors, generally referred to as geophones when used on land, and as hydrophones when used offshore. On land the source of seismic energy can be a high explosive charge electrically detonated in a borehole located at a selected point on a terrain, or another energy source having capacity for delivering a series of impacts or mechanical vibrations to the earth""s surface. The acoustic waves generated in the earth by these sources are transmitted back from strata boundaries and reach the surface of the earth at varying intervals of time, depending on the distance and the characteristics of the subsurface traversed. These returning waves are detected by the geophones, which function to transduce such acoustic waves into representative electrical signals. In use an array of geophones is generally laid out along a line to form a series of observation stations within a desired locality, the source injects acoustic signals into the earth, and the detected signals are recorded for later processing using digital computers where the data are generally quantized as digital sample points such that each sample point may be operated on individually. Accordingly, seismic field records are reduced to vertical and horizontal cross sections which approximate subsurface features. The geophone array is then moved along the line to a new position and the process repeated to provide a seismic survey. More recently seismic surveys involve geophones and sources laid out in generally rectangular grids covering an area of interest so as to expand areal coverage and enable construction of three dimensional (3D) views of reflector positions over wide areas.
The general principle of offshore prospecting consists of using a seismic source to create a disturbance in the marine environment (e.g. by releasing a volume of air or steam into the water, by varying the volume of an immersed body, by an implosion, etc), and in using hydrophone detectors towed by a prospecting ship, or by geophones placed on the seabed, to obtain seismic data for extracting useful information concerning the geology of the subsoil. Typically in offshore exploration, a long streamer cable that electrically connects multiple hydrophones is towed behind the prospecting ship. Acoustic waves are periodically generated in the water, reflected from subterranean earth layers, and detected by the hydrophones. The hydrophones convert the detected seismic waves into representative electrical signals, which may be processed on the ship, and/or recorded on a storage medium such as magnetic tape for later processing.
As oil and gas production from known reservoirs and producing provinces declines, explorationists seek new areas in which to find hydrocarbons. Many of the new areas under exploration contain complex geological structures that are difficult to image with 2D techniques. Accordingly, 3D seismic processing has come into common use for mapping subterranean structures associated with oil and gas accumulations. Geophysicists, however, are well aware that a 2D seismic record section or 3D view is not a true reflectivity from the earth, but is instead a transformation of the earth""s reflectivity into a plane where each recorded event is located vertically beneath the source/receiver midpoint. Steep dip, curved surfaces, buried foci, faults and other discontinuities in subterranean structure each contribute their unique characteristics to the seismic record and, in complexly faulted and folded areas, make interpretation of the geological layering from the seismic record extremely difficult. Migration is the inverse transformation that carries the plane of recorded events into a true 3D reflectivity of the earth, thereby placing reflections from dipping beds in their correct location and collapsing diffractions.
Of the various available migration methods, wave-equation migration is considered to be superior because it is based on accurate propagation of seismic waves through complex models of the earth. Wave equation computations using numerical techniques have led to a procedure called reverse time migration (RTM). By this procedure, the wavefield recorded at the surface is imaged in depth using a model of earth velocities in a numerical solution of the wave-equation. The wavefield at the surface is used as a boundary condition for the numerical computations. Proceeding by inserting the data at the surface of a computational grid for each record time step, starting with the last recorded time sample, and ending with the first, the wavefield migrates to the position from which the reflections originated.
Recently, numerical finite-difference methods have been widely used for obtaining seismic images based on pre-stack and post-stack RTM techniques, and for obtaining synthetic seismograms for studying exploration related problems. A practical procedure for doing reverse time migration is disclosed in a publication, Mufti, I. R., et al, xe2x80x9cFinite-Difference Depth Migration of Exploration-Scale 3D Seismic Data,xe2x80x9d Geophysics, Vol. 61. No. 3 (May-June 1996), which is incorporated herein by reference. RTM, which requires enormous computer resources as compared to simpler or less accurate migration algorithms, has recently been applied to 3D seismic data. An improved image results from accuracy (dynamic as well as kinematic) of the finite-difference method over conventional normal-moveout and raytracing based seismic imaging methods. In this RTM procedure a finite-difference earth model, which is based on the best estimate of subsurface velocities, is required. This involves dividing the model simulation space into a large number of elementary grid blocks, and assigning a velocity value to each grid.
A common problem with the finite-difference migration method is that simulation of wave propagation in an extensive portion of the earth must be modeled with limited computer resources, i.e., mainly limited central RAM memory in the computer such that the computational size of the model is truncated due to limited computer memory size. Accordingly, artificial boundaries that act as perfect reflectors are imposed by computer memory limitations, and if not properly handled cause unwanted reflection artifacts. These unwanted reflection artifacts may disrupt the image in complex ways and are not easily removed after the modeling and imaging processes are carried out.
One method of absorbing advancing waves in a 3D computational grid is an extensively used method known as a one-way wave equation (1WWE), proposed by Clayton and Engquist, 1980 xe2x80x9cAbsorbing Boundary Conditions for Wave-Equation Migrationxe2x80x9d, Geophysics, Vol. 45 pp 895-904, and Keys, 1985 xe2x80x9cAbsorbing Boundary Conditions for Acoustic Mediaxe2x80x9d, Geophysics, Vol. 50 pp 89-902. This 1WWE method, however, is effective for absorbing advancing waves arriving only at or near a selective incident angle. More recently, much attention has been given to a method called the perfectly matched layer (PML) method, which was developed by Berenger (1994) xe2x80x9cA Perfectly Matched Layer for the Absorpiton of Electromatic Wavesxe2x80x9d, Journal of Computational Physics Vol. 114, pp 185-200, for use in a finite-difference time domain computations. The PML method, which was first applied to a 2D transverse electric wave problem, can absorb waves propagating at both near normal and high incident angles, and thus provides a reflectionless interface between the propagation region of interest and the PML layers at all incident angles of arriving waves. A PML layer can be conveniently applied on a planar face. However an edge or a corner of a rectangular finite difference numerical model volume requires a special tensor, which adds significantly to computational requirements. Yuan, et al., xe2x80x9cFormulations and Validation of Berenger""s PML Absorbing boundary for the FDTD Simulation of Acoustic Scattering,xe2x80x9d IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 44, pp 816-822, July 1997, lays out a detailed implementation of the PML medium suitable for both 2D and 3D acoustic applications. However, this implementation of the PML medium requires computations for at least six component variables, and accordingly is an expensive approach in terms of computer capacity when applied to large 3D seismic problems.
Accordingly, an urgent need exists for a cost effective and straightforward method to reduce boundary reflection artifacts in finite-difference 3D wave equation computations.
The basic computing requirements for successful 3D finite-difference acoustic wave equation seismic migration computations include a large memory and a high-speed processor, or preferably a multiprocessor massively parallel computing system for execution of numerical simulation software.
It is an object of the present invention to provide an absorbing boundary region for termination of a finite-difference time domain computation region.
A more specific object of the present invention is to eliminate unwanted artifacts (i.e. reflections and diffractions) arising from truncated boundaries imposed on a numerical acoustic model.
A still more specific object is to absorb simulated wavefield energy arriving at arbitrary incident angles at the boundaries of a finite-difference computational model, without significantly increasing computer memory capacity and/or computation time.
A further object is to reduce computational time and computer memory requirements by reducing the number of variables and computational layers that are employed in the absorbing boundaries.
A still further object of this invention is to produce a computer program which generates high resolution images of seismic wave propagation.
According to the present invention the foregoing and other objectives and advantages are attained in a method of eliminating reflection artifacts which are introduced in finite-difference acoustic wavefield propagation computations by reflections from artificial boundaries. These artificial boundaries are imposed by a memory size lacking capacity for containing an extensive earth velocity model based on a large 3D volumetric seismic survey. In this method, two types of absorbing boundaries conditions are imposed to form an absorbing region surrounding the truncated earth velocity model boundaries, thus providing an absorbing region which presents an essentially reflectionless interface for outgoing seismic waves. Accordingly, the truncated boundaries are terminated in a manner that absorbs simulated wave energy contained in an advancing seismic wave, and unwanted reflection artifacts, which would otherwise be reflected from artificial boundaries, are eliminated.
The reflection elimination method according to this invention uses a hybrid modeling scheme to provide the absorbing region which is effective at essentially all incidence angles of the advancing wave. PML conditions, which can absorb waves at essentially all incidence angles, are coupled to finite-difference wavefield propagation computation regions at vertical side faces and the bottom of rectangular shape memory volumes. In cooperation with the imposed PML conditions, 1WWE conditions, which absorb advancing waves at or near a selective incidence angle, are coupled to the exterior surface of the PML absorbing layers, and to the edges and corners of the truncated earth velocity model, where implementing the PML type conditions is significantly more complicated when compared to the flat surface PML implementation.
When using the PML and 1WWE conditions in combination with the finite-difference computations, the procedure for computing wave propagation calls for updating variables in the various regions in the following order: first the time and space variables in the finite-difference region, followed by applying 1WWE conditions at the vertical edges and corners, updating particle velocity and a specifically defined variable, q, in the interior of the PML region, applying 1WWE at the exterior of the PML regional, updating pressure in the interior of the PML, and updating pressure in the two layers in the wavefield propagation region immediately above the absorbing region.
In another aspect of this invention a highly efficient representation of the PML absorbing conditions is achieved by reducing the number of variable components used in the PML computations from a conventional six components to four components. This reduction is achieved based on introducing a new variable, herein identified as xe2x80x9cqxe2x80x9d, which is defined in a manner that eliminates two of the usual six variables, as will be illustrated hereinafter.
In yet another aspect of the invention, apparatus comprises a computer programmed for migration of seismic data which eliminates undesired reflections and diffractions from truncated boundaries imposed on 3D finite-difference acoustic applications. For a large scale 3D survey, the computer is a massively parallel machine.
The hybrid approach for implementing absorbing regions is beneficial for its rapid attenuation of waves arriving at arbitrary angles, which leads to elimination of reflection artifacts, and is ideal for implementation on massively parallel super computers such as a CRAY T3E or an IBM SP2. In these large 3D applications, where the computational regions for the velocity model are desirably partitioned between the multiple processors, artificial boundaries are inherently created between the multiple processors. Accordingly, employing the hybrid modeling scheme of this invention optimizes the absorbing region to reduce memory usage and the related computational costs, especially when applied to large 3D seismic surveys.